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We derive envelope equations which generate the Maxwell-Lorentz model and 
describe the interaction of optical pulses with plasmonic oscillations in metal 
nanoparticle composites. A family of solitary wave solutions is found which 
is analogous to self-induced transparency in Maxwell-Bloch. The evolution of 
incident optical pulses is studied numerically as are the collision dynamics of 
the solitary waves. These simulations reveal that the collision dynamics vary 
from near perfectly elastic to highly radiative depending on the relative phase 
of the initial pulses, (c) 2008 Optical Society of America 

OCIS codes: 190.4400, 250.5530, 260.3910 

Quantum effects in metal nanoparticles driven by a resonant optical field play an impor- 
tant role in inducing a strong nonlinear response, as was recently shown. ^'^ In this Letter we 
consider the nonlinear resonant interaction of ultrashort optical pulses with metal nanopar- 
ticles distributed uniformly in a host medium. We restrict to the case of composite materials 
for which the resonance frequencies of the host medium are well separated from those of the 
nanoparticles. Examples include silver or gold spherical or spheroidal nanoparticles embed- 
ded in Si02. In these cases, the plasmonic resonance frequencies are in the visible part of the 
spectrum while the resonance of the host is in the ultraviolet. 

Light interaction with metal nanoparticles can be described by a system consisting of 
Maxwell's equations for the electric field, and an oscillator equation describing the displace- 
ment of conduction electrons in the metal nanoparticles from equilibrium (plasmonic oscil- 
lations). The nanoparticles are much smaller than the optical carrier wavelength Aq. This 
allows light scattering and spatial effects in the nanoparticles to be neglected. As shown by 
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Rautian^ and Drachev, et al.,^ who further developed the earher work by Hache, et aL,^ 
the response of the conduction electrons in the metal nanoparticles to an external electric 
field induces a leading-order cubic nonlinearity. The interaction of the electric field with 
plasmonic oscillations in nanoparticles with resonance frequency uj^ in the presence of this 
cubic nonlinearity can be described by the forced Duffing equation 

Qtt + ojIQ + xQ^ = {e/m)£. (1) 

In this expression Q represents plasmon displacement from equilibrium, T is time, x is the 
coefficient of nonlinearity, e and m are the electron charge and rest mass, respectively, and 
S is the electric field. The tilde is used to denote rapidly-varying quantities. The nonlinear 
coefficient x can be estimated by comparing the off-resonance nonlinear response in Eq. (Q) 
with the Drude nonlinearity for (non-resonant) conduction electrons in metal nanoparticles. 
That susceptibility is characterized by^ x^^^ — Ne^a? / {mh^ujQ). Here a and are the radius 
of the nanoparticle and the conduction electron density of the metal, respectively, and ujq is 
the optical carrier frequency. This results in the estimate k ~ {maujQ/hY. 

We are interested in pulse dynamics which vary on a much slower scale than the plasmonic, 
host atom, and carrier wave oscillations, and can be described using a slowly- varying envelope 
approximation. In this approximation, Eq. becomes 

^Qt + (c^r.-c^o)Q+(3x/2cuo) I Q = -(e/2ma;o)^, (2) 

where the slowly-varying envelopes of the electric field and plasmonic oscillations are rep- 
resented by £ and Q respectively. Maxwell's equation couples to the material polarization 
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induced by the plasmonic oscillations. The equation for the electric field envelope is 



where Z is the propagation coordinate, Vg is group velocity, c is the speed of light, no is the 
refractive index evaluated at the carrier frequency c^o? is the product of the conduction 
electron density and the metal filling factor p (the fraction of the composite occupied by 
metal). Na is the concentration of host atoms, d is the projection of the dipole matrix element 
in the direction of the electric field polarization, and = — cuq is detuning from the 
resonance frequency of host atoms. The last two terms in Eq. (0) represent corrections to the 
refractive index and group index due to the off-resonance interaction with the host medium, 
which, for illustration, is considered as an ensemble of two-level atoms. This equation is 
derived from the Maxwell-Bloch equations in the non-resonant case by considering as a 
large parameter and applying the adiabatic following approximation.^ Additional resonances 
would produce similar terms. We consider the case where optical pulse intensity and duration 
as well as composite material parameters are such that the characteristic length of resonant 
light interaction with plasmonic oscillations is much smaller than the characteristic lengths 
for both group velocity dispersion and nonlinearity induced by the host medium. Therefore 
the terms representing these effects are omitted from Eq. 

In a composite material, the sizes and shapes of metal nanoparticles vary due to limited 
fabrication tolerances. It is known that the plasmon resonance in spherical metal nanopar- 
ticles depends weakly on size in the range between 10 and 50nm,^ so that variations in size 
are not important. However, variations in the shape and orientation of the nanoparticles 
can significantly change plasmonic resonance frequencies. This results in a broadening of the 
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resonance line of the bulk composite. The angle brackets {Q(t, z, u)) = Q{t, z, uj)g{uj)duj 
denote averaging over the distribution g{uj) of the resonance frequencies (line shape). Defin- 
ing 

£ = -E ^-J exp{iKZ), Q = Quo\ —exp{iksZ), (4) 

e V ox V ox 

where ks = 27rujQNa\d\'^ / criohAa , and introducing the copropagating coordinate system z = 

{ujp/4:cnQLjQ)Z, t = uJo(T — Z/u), {u here is shifted group velocity defined as = v^^ + 

{27rLJoNa\d\'^ / cnohAl) , oj^ = AirNpe^/m, u = {ur — ujq)/ujq, Eqs. (0) and Q can be reduced 

to the simpler form 

iE, = {Q), iQt + ujQ + \Q\''Q = E. (5) 

These equations represent a generalization of the classical Maxwell-Lorentz model. In the 
case of identical nanoparticles, the averaging in (0) can be reduced to a single dimensionless 
frequency u [i.e. detuning frequency distribution g{uj) = S{uj — cu)]. Under this condition the 
system has solitary wave solutions: 

p., X v'^^ exp [i^ + int - iKj - ixiC)] exp|-2zx(0] 



where i = [z - v{t - r)]/^, x(0 = arctan [F tanh (e/2eo)], = 1/2(1 - K-^f'^ T = 
[(1 — K)/{1 + i^)]^/^, and K = {uj — VL)/2^/v. These solutions are parameterized by velocity 
w, frequency 1), phase shift y?, and position r. The velocity v is the amount by which the 
wave is slowed from the copropagating frame velocity u. Thus in the laboratory frame, the 
actual pulse velocity is u — v. The quantity .^o must be real, hence 1 — > 0. Thus the 
condition for existence of these solutions is la) — fij < 2^/v. This choice of parameters provides 
relatively simple mathematical expressions for the solitary waves. In practice it is easier to 



both control and measure peak amplitude, A = 2f^/^(l — Ky^'^, than the pulse velocity, 
therefore A, Q, (f, r form a more suitable set of parameters. Given the pulse amplitude 
A, the corresponding velocity parameter depends on the value of the quantity oj — VL. If 
uj = VL^ then v = (74/2)^/'^ trivially. For the case when co ^ il, write the amplitude as 
A = 2v^^^{l — a\iu — ^l\/2^/vY^'^, where the parameter a = sgn(c(; — fi). Then defining 
V = (2y^/\iu — il\ —a)^^"^ and A = 27 /2\0 — A leads to an expression for the rescaled 
velocity v = {y~^^^ — cry^^^) /^/S, where y = o"[(A^ + cr)^/^ — A]^/^. In this calculation, the 
appropriate branches have been chosen so that the expressions are consistent with reality 
and positivity conditions on the parameters. 

In optics it has become standard practice to refer to certain solutions of nonintegrable 
systems as solitons. These solutions are characterized as solitary waves which are robust to 
external perturbations including collisions with other solitary waves. In addition, arbitrary 
initial data for these "soliton" supporting systems tends to evolve into a sum of solitary waves 
and continuous radiation. The remainder of this Letter details our numerical investigation of 
these properties, in which Eqs. (0) are integrated in the case of delta-distributed resonance 
frequencies and zero detuning. 

Numerical simulations of the evolution of Gaussian initial data are presented in Figs. ^ 
and b. In Fig.[T^ the initial condition E(t, 0) = exp(— 1^/2)/2 simply evolves into continuous 
radiation, while in Fig. QJd the initial condition E{t,0) = 2exp(— 1^/2) emits some radiation 
but also achieves energy confinement and persists as a soliton. This behavior is similar to the 
self-induced transparency exhibited by the Maxwell-Bloch equations^ which describe optical 
pulse interaction with resonant two-level media. As the amplitude of the initial pulse is 
increased, the pulse splits into two [illustrated in Figl^ti- for the initial condition E{t,0) = 
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5exp(— or more solitons and emits continuous radiation. The inset shows a comparison 
between the numerics and the analytic form of the sohtary waves of Eq. The velocity 
and frequency of the solitary waves are obtained from measurements of their amplitudes 
and half widths in our simulation. The agreement of the analytic forms with the simulation 
results indicates that the system self-selects the solitary waves presented in Eq. Fig. |2b 
shows the amplitude (s) of these solitons as a function of input pulse amplitude. As the input 
pulse amplitude increases, the output soliton amplitude also increases until a bifurcation 
occurs and a new soliton emerges. Increasing the input pulse amplitude further results in 
the production of more solitons, along with continuous radiation. 

Two simulations illustrating collision dynamics are presented in Figs.Eti' andEb, where the 
sum of two well-separated solitary waves is used as the initial condition. Fig. illustrates 
an in-phase collision, in which the relative phase Aip = cpi — (p2 = 0, where the subscripts 
identify the soliton. The other initial soliton parameters are fi = 1, fii = 0, ri = —10, 
V2 = 2, Q2 = 0, and T2 = 0. In Fig. an out-of-phase A(p = it collision is illustrated. The 
same parameters are used except for a shift in the relative phase. This results in a much 
faster collision. In both simulations the solitary waves persist after the interaction, although 
their characteristic parameters undergo shifts and radiation is emitted during the collision. 
A numerical study indicates that the collisions are quasielastic for values in the approximate 
interval A(/9 G (vr, 27r). For some A(/9 values away from this interval, simulations show that 
one of the solitons is completely destroyed while the other persists. A detailed analysis of 
the dependence on initial parameters will be presented elsewhere. 

The soliton phenomena described above occurs at light intensities such that the dimen- 
sionless field amplitude E is at least of order one. The intensity at which E is order one 
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estimated as / ~ (c/x)(mc(jQ/e)^ ^ c{hjjQ/ eaY . For the particle radius a = 20nm and the 
carrier wavelength Aq = 500nm, this results in light intensity of / ~ lOGW/cm^ which can be 
easily obtained with ultrashort laser pulses. The optical pulse durations for which this model 
is valid are limited by the condition Au <^ uq, (the spectral width of the pulse must be much 
smaller than the carrier frequency) required by the slowly varying envelope approximation. 
The pulse duration should also be much shorter than the characteristic plasmonic oscillation 
damping time, which is determined by the time required for electron thermalization in the 
metal nanoparticles (~400fs^'^). The envelope approximation is appropriate for pulses with 
width T > 20fs. 

In summary, a family of solitary wave solutions is derived in the envelope approximation for 
the Maxwell wave and Duffing oscillator equations, showing that energy confinement is pos- 
sible for resonant optical pulse interaction with plasmonic oscillations in metal nanoparticles. 
The existence condition for these solutions is presented. Numerical simulations show that 
stable solitary waves evolve from arbitrarily-shaped initial pulses with sufficient amplitudes 
and exhibit behavior analogous to self-induced transparency in Maxwell- Bio ch. Simulations 
also reveal that the collision dynamics are highly dependent on initial soliton parameters, 
behaving quasielastically in some regimes but having radically different behavior in others. 
The authors are grateful to V. P. Drachev for helpful discussions. In addition, we would like 
to acknowledge funding under Arizona Proposition 301, LANL, and NSF. 
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List of Figure Captions 

n Evolution of electric field amplitude with initial conditions exp(— 1^/2)/2 (left) and 
2exp(-tV2) (right). 

12] Left: Evolution of electric field amplitude with initial condition 5exp(— 1^/2). The inset 
shows a comparison of the numerics (dashed line) with the analytic form of the solitary wave 
solutions (solid line). Right: output solitary wave amplitude(s) as a function of Gaussian 
input pulse amplitude Aq, where the initial condition is given by exp(— 1^/2). 
OJElectric field amplitude showing collision dynamics of solitons for different values of relative 
phase. Left: Ay9 = 0; right: Ay^ = tt. 
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Fig. 1. Evolution of electric field amplitude with initial conditions 
exp(-tV2)/2 (left) and 2exp(-tV2) (right). 
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Fig. 2. Left: Evolution of electric field amplitude with initial condition 
5exp(— The inset shows a comparison of the numerics (dashed line) 
with the analytic form of the solitary wave solutions (solid line). Right: out- 
put solitary wave amplitude(s) as a function of Gaussian input pulse amplitude 
Ao, where the initial condition is given by exp(— 1^/2). 
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